library(foreign)
library(polycor)

setwd("~/Google Drive/Broockman/Completed/Extremity vs Consistency/ReplicationData/PanelSurveyApril2014")

all.data <- read.dta("ssi-spatial-all-merged-with-wave1.dta")
names(all.data)

####Figure
make.test.retest.graph <- function(data, filename){
	save <- TRUE
	names <- c("Environment","Social\nSecurity","Guns","Health","Immigration","Taxes","Abortion","Medicare","Gay Rights","Unions","Contraception","Education")
	cat(cbind(names(data)[1:12],names)) #spot check
	cat(cbind(names(data)[13:24],names)) #spot check
	
	col.ramp <- colorRampPalette(c("red", "green"))
	
	start.index.pre <- 0
	start.index.post <- 12
	scat <- function(pre.issue, post.issue){
		x <- data[,start.index.pre + pre.issue]
		y <- data[,start.index.post + post.issue]
		corr <- polychor(x,y)
		corr.rounded <- round(corr*100)
		if(corr.rounded < 1) corr.rounded <- 1 #for color coding figures
		if(corr.rounded > 70) corr.rounded <- 70
		plot(x + rnorm(nrow(data), sd=.1), y + rnorm(nrow(data), sd=.1),
			xlab = "",#paste0("First Wave: ", names[pre.issue]),
			ylab = "",#paste0("Second Wave: ", names[post.issue]),
			pch = 16,
			cex = .3,
			axes=FALSE
		)
		if(pre.issue == post.issue) box(which="figure")
		title(sprintf("%.2f",corr),
			col.main = col.ramp(70)[corr.rounded]
			)
		lines(loess.smooth(x,y), lwd=4, col='red')
		return(corr)
	}
	
	lab <- function(index, pre){
		if(pre) lab <- paste0(names[index], "\nFirst Test")
		if(!pre) lab <- paste0(names[index], "\nRetest")
		plot(-1,-1,xlim=c(0,1),ylim=c(0,1),axes=FALSE,xlab="",ylab="")
		text(0.5,0.5, lab, cex=.7, srt=45,font=2)
	}
	
	issues <- 1:12
	n.issues <- length(issues)
	corr.mat <- matrix(nrow=max(issues),ncol=max(issues))
	
	if(save) png(filename, width=2500, height=2500, res=330)
	par(mfrow=c(n.issues+1,n.issues+1), mar = c(0,0,1,0))
	for(i in c(0,issues)){
		for(j in c(0,issues)){
			if(i > 0 & j > 0) corr.mat[i,j] <- scat(j,i)
			if(i == 0 & j > 0) lab(j, TRUE)
			if(j == 0 & i > 0) lab(i, FALSE)
			if(i == 0 & j == 0) plot.new()
		}
	}
	if(save) dev.off()
	
	#####Average correlations
	corr.mat.no.diag <- corr.mat
	for(i in 1:nrow(corr.mat)) corr.mat.no.diag[i,i] <- NA
	cat(mean(corr.mat.no.diag, na.rm=TRUE))
	
	corrs <- vector()
	for(i in 1:nrow(corr.mat)) corrs <- c(corrs,corr.mat[i,i])
	cat(mean(corrs))
}

pkscale <- all.data$pkscale
multis <- all.data[,11+c(1:12,14:25)]
make.test.retest.graph(multis, "testretest-graph.png")

above.pk.median <- multis[pkscale>median(pkscale),]
below.pk.median <- multis[pkscale<median(pkscale),]
make.test.retest.graph(above.pk.median, "testretest-graph-highknowledge.png")
make.test.retest.graph(below.pk.median, "testretest-graph-lowknowledge.png")